'''Analyze Quintara Sequencing output .fas files
	Conventions: Samples are S#.N, where # is a construct, N is a clone
	Inputs: Construct Map, directory of sequencing files, (optional) thresholds for good reads.
	Output: Directory of high signal reads for simultaneous alignment in ApE
	WishList: Output: MUSCLE alignments of reads, +/- per clone, explanation of issues.'''

from Bio import Seq,SeqRecord,SeqIO,Align,AlignIO
from Bio.Align.Applications import MuscleCommandline
import csv
import string
import datetime
import os
import pdb
		
def today():
    '''YYYY-MM-DD time to MMDDYY'''	
    date_str = str(datetime.datetime.today())
    return date_str[5:7]+date_str[8:10]+date_str[2:4]

def sequencing(seq_dir='../Dropbox/Thomas Rotation/Sequencing/2012-08-19 pSNV/',lower=50,upper=850,minlen=500):
    '''IO Wrapper'''
    ##########
    #pdb.set_trace()
    ##########

    def parse_construct():
        '''Construct to new FASTA Collection to align'''
        
        def add_samples():
            '''Add Sequencing Samples.
                Outputs: path to collection file'''
            seq_list_dir = seq_dir+in_id+'/'
            if not os.path.exists(seq_list_dir):
                os.mkdir(seq_list_dir)
            for filename in os.listdir(seq_dir):
                if filename.find('fas')>-1:
                    #Look in all quintara .fasta files
                    for seq_record in SeqIO.parse(seq_dir+filename,'fasta'):
                        #Check for current sample_id (case-insensitive)
                        if seq_record.id.upper().find(sample_id.upper())>-1:
                            #Read must exceed minimum length
                            if len(seq_record.seq)>minlen:
                                #Add only lower:upper of that string, unless upper exceeds read length
                                seq_list.append(seq_record[lower-1:min(upper,len(seq_record.seq))])
                    ##########
                    #pdb.set_trace()
                    ##########
                    for i,seq_record in enumerate(seq_list):
                        if i==0:
                            #keep construct copy in genbank format for ApE display
                            SeqIO.write(seq_record,seq_list_dir+seq_record.id+'.gb','genbank')
                        else:
                            SeqIO.write(seq_record,seq_list_dir+seq_record.id+'.fasta','fasta')
            return seq_list
        ##########
        #pdb.set_trace()
        ##########		
        in_id = sample_path.split('/')[-1].split('.')[0]
        seq_construct = SeqIO.read(sample_path, "genbank")
        seq_construct.id = 'TEMPCOPY '+in_id
        seq_list.append(seq_construct)
        return add_samples()


    key = csv.reader(open(seq_dir+'key.csv','rU'), delimiter=',')
    for i,row in enumerate(key):
        if i>0 and len(row[2])>1:
            sample_id = row[0]
            sample_path = row[2]
            seq_list = []
            parse_construct()
        
        ''' MUSCLE SECTION
        if out_file:
            muscle_cline = MuscleCommandline(input=in_file,out=out_file)
        else: 
            muscle_cline = MuscleCommandline(input=in_file)
        stdout, stderr = muscle_cline()
        '''

def main(*argv):
    return 'main'
	
if __name__=='__main__':
    main()


